---
title: "7-2 WHO 8-point Outcome Scale"
description: |
Analyses of the WHO 8-point outcome scale outcome.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
toc-depth: 5
---
```{r}
#| label: pkgs
#| code-summary: Load packages
library (ASCOTr)
library (tidyverse)
library (lubridate)
library (kableExtra)
library (patchwork)
library (cmdstanr)
library (posterior)
library (bayesplot)
library (ggdist)
library (lme4)
library (broom)
library (broom.mixed)
library (bayestestR)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
all_data <- readRDS (file.path (ASCOT_DATA, "all_data.rds" ))
# Create censored D28 outcome for the two participants
all_data <- all_data |> mutate (
D28_who_low = case_when (
StudyPatientID %in% c ("ALF00006" , "ALF00012" ) ~ 1L,
TRUE ~ D28_who
),
D28_who_high = case_when (
StudyPatientID %in% c ("ALF00006" , "ALF00012" ) ~ 2L,
TRUE ~ D28_who
)
)
# FAS-ITT
fas_itt_dat <- ASCOTr::: make_fas_itt_set (all_data)
fas_itt_nona_dat <- fas_itt_dat |>
filter (! is.na (D28_who))
# ACS-ITT
acs_itt_dat <- ASCOTr::: make_acs_itt_set (all_data)
acs_itt_nona_dat <- acs_itt_dat |>
filter (! is.na (D28_who))
# AVS-ITT
avs_itt_dat <- ASCOTr::: make_avs_itt_set (all_data)
avs_itt_nona_dat <- avs_itt_dat |>
filter (! is.na (D28_who))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative" ), dir = "stan" )
ordmodcens0 <- compile_cmdstanr_mod (
file.path ("ordinal" , "logit_cumulative_cens" ), dir = "stan" )
ordmod_site <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_site" ), dir = "stan" )
ordmod <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_epoch_site" ), dir = "stan" )
ordmodcens <- compile_cmdstanr_mod (
file.path ("ordinal" , "logit_cumulative_epoch_site_cens" ), dir = "stan" )
logistic <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site_epoch" ), dir = "stan" )
```
```{r}
#| label: helpers
#| code-summary: Helper functions
make_primary_model_data2 <- function (dat,
outcome = "PO" ,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_int = 2.5 ,
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
beta_sd_trt = 1 ,
ctr = contr.equalprior,
includeA = TRUE ,
includeC = TRUE ,
intercept = TRUE ,
...) {
X <- make_X_design (dat, vars = vars, ctr = ctr, includeA = includeA, includeC = includeC, intercept = intercept)
nXtrt <- sum (grepl ("rand" , colnames (X)))
epoch <- dat[["epoch" ]]
M_epoch <- max (epoch)
region <- dat[["ctry_num" ]]
M_region <- max (region)
site <- dat[["site_num" ]]
M_site <- max (site)
region_by_site <- region_by_site <- dat |>
dplyr:: count (ctry_num, site_num) |>
pull (ctry_num)
y_raw <- as.matrix (dat[, outcome])
y_mod <- y_raw
N <- dim (X)[1 ]
K <- dim (X)[2 ]
beta_sd <- c (rep (beta_sd_trt, nXtrt), beta_sd_var)
if (intercept) {
beta_sd <- c (beta_sd_int, beta_sd)
}
out <- list (
N = N, K = K, X = X, y = y_mod, y_raw = y_raw,
J = max (y_mod), p_par = rep (2 / max (y_mod), max (y_mod)),
M_region = M_region, region = region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
beta_sd = beta_sd
)
if (includeA) {
out <- c (out, list (ctrA = attr (X, "contrasts" )[["randA" ]]))
}
if (includeC) {
out <- c (out, list (ctrC = attr (X, "contrasts" )[["randC" ]]))
}
return (out)
}
fit_primary_model2 <- function (dat = NULL ,
model = NULL ,
outcome = "PO" ,
vars = c ("inelgc3" , "agegte60" , "sexF" , "supp_oxy2" , "ctry" ),
beta_sd_int = 2.5 ,
beta_sd_var = c (10 , 2.5 , 2.5 , 2.5 , 1 , 1 ),
beta_sd_trt = 1 ,
ctr = contr.equalprior,
includeA = TRUE ,
includeC = TRUE ,
intercept = TRUE ,
seed = 32915 ,
adapt_delta = 0.99 ,
iter_sampling = 2500 ,
chains = 8 ,
...) {
mdat <- make_primary_model_data2 (
dat,
outcome = outcome,
vars = vars,
beta_sd_int = beta_sd_int,
beta_sd_var = beta_sd_var,
beta_sd_trt = beta_sd_trt,
ctr = ctr,
includeA = includeA,
includeC = includeC,
intercept = intercept,
...
)
snk <- capture.output (
mfit <- model[["sample" ]](
data = mdat,
seed = seed,
adapt_delta = adapt_delta,
refresh = 0 ,
iter_warmup = 1000 ,
iter_sampling = iter_sampling,
chains = chains,
parallel_chains = min (chains, parallel:: detectCores ())
)
)
mpars <- mfit$ metadata ()$ model_params
keep <- mpars[! grepl ("(_raw|epsilon_)" , mpars)]
mdrws <- as_draws_rvars (mfit$ draws (keep))
names (mdrws$ beta) <- colnames (mdat$ X)
if (any (grepl ("site" , keep))) {
site_map <- dat %>% dplyr:: count (site_num, site)
names (mdrws$ gamma_site) <- site_map$ site
}
if (any (grepl ("epoch" , keep))) {
epoch_map <- dat %>% dplyr:: count (epoch, epoch_lab)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
}
if (includeA) {
Ca <- mdat$ ctrA
mdrws$ Acon <- as.vector (Ca %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))])
names (mdrws$ Acon) <- rownames (Ca)
mdrws$ Atrt <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ AOR <- exp (mdrws$ Atrt)
}
if (includeC) {
Cc <- mdat$ ctrC
mdrws$ Ccon <- as.vector (Cc %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))])
names (mdrws$ Ccon) <- rownames (Cc)
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ COR <- exp (mdrws$ Ctrt)
}
mdrws$ OR <- exp (mdrws$ beta[! grepl ("(Intercept|rand)" , names (mdrws$ beta))])
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
make_summary_table_anticoagulation <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (
CAssignment = factor (CAssignment, labels = intervention_labels2 ()$ CAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (D28_who)),
Deaths = sprintf ("%i (%.0f%%)" ,
sum (D28_who == 8 , na.rm = TRUE ),
100 * mean (D28_who == 8 , na.rm = TRUE )),
` Hospitalised ` = sprintf ("%i (%.0f%%)" ,
sum (D28_who > 2 & D28_who < 8 , na.rm = TRUE ),
100 * mean (D28_who > 2 & D28_who < 8 , na.rm = TRUE )),
` WHO, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (D28_who, na.rm = T),
quantile (D28_who, 0.25 , na.rm = TRUE ),
quantile (D28_who, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (D28_who)),
Deaths = sprintf ("%i (%.0f%%)" ,
sum (D28_who == 8 , na.rm = TRUE ),
100 * mean (D28_who == 8 , na.rm = TRUE )),
` Hospitalised ` = sprintf ("%i (%.0f%%)" ,
sum (D28_who > 2 & D28_who < 8 , na.rm = TRUE ),
100 * mean (D28_who > 2 & D28_who < 8 , na.rm = TRUE )),
` WHO, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (D28_who, na.rm = T),
quantile (D28_who, 0.25 , na.rm = TRUE ),
quantile (D28_who, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_summary_table_antiviral <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (
AAssignment = factor (AAssignment, labels = intervention_labels2 ()$ AAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (D28_who)),
Deaths = sprintf ("%i (%.0f%%)" ,
sum (D28_who == 8 , na.rm = TRUE ),
100 * mean (D28_who == 8 , na.rm = TRUE )),
` Hospitalised ` = sprintf ("%i (%.0f%%)" ,
sum (D28_who > 2 & D28_who < 8 , na.rm = TRUE ),
100 * mean (D28_who > 2 & D28_who < 8 , na.rm = TRUE )),
` WHO, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (D28_who, na.rm = T),
quantile (D28_who, 0.25 , na.rm = TRUE ),
quantile (D28_who, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (AAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (D28_who)),
Deaths = sprintf ("%i (%.0f%%)" ,
sum (D28_who == 8 , na.rm = TRUE ),
100 * mean (D28_who == 8 , na.rm = TRUE )),
` Hospitalised ` = sprintf ("%i (%.0f%%)" ,
sum (D28_who > 2 & D28_who < 8 , na.rm = TRUE ),
100 * mean (D28_who > 2 & D28_who < 8 , na.rm = TRUE )),
` WHO, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (D28_who, na.rm = T),
quantile (D28_who, 0.25 , na.rm = TRUE ),
quantile (D28_who, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Antiviral \n intervention ` = AAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
```
# Descriptive {#descriptive}
## Anticoagulation
```{r}
#| label: tbl-7-2-summary-anticoagulation
#| code-summary: Summary of WHO outcome by arm
#| tbl-cap: Summary of WHO scale at day 28 by treatment group, anticoagulation domain, FAS-ITT.
save_tex_table (
make_summary_table_anticoagulation (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-2-anticoagulation-summary" ))
make_summary_table_anticoagulation (fas_itt_dat)
```
```{r}
#| label: fig-who-anticoagulation-fas-itt
#| fig-cap: Day 28 WHO score, anticoagulation domain, FAS-ITT.
p1 <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
group_by (CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (CAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Anticoagulation intervention" )
p2 <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (CAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-2-anticoagulation.pdf" ), p2, height = 3 , width = 6 )
p2
```
```{r}
#| label: fig-who-anticoagulation-avs-itt
#| fig-cap: Day 28 WHO score, anticoagulation domain, AVS-ITT.
p2 <- avs_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (CAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-2-anticoagulation-avs-itt.pdf" ), p2, height = 3 , width = 6 )
p2
```
## Antiviral
There were two participants randomised to the antiviral domain who had unknown day 28 status on the WHO scale.
For these two participants it was known that they were alive and out of hospital, but not whether they were grade 1 or grade 2.
In the models, these participants are treated as censored.
```{r}
#| label: tbl-7-2-summary-antiviral
#| code-summary: Summary of WHO outcome by arm
#| tbl-cap: Summary of WHO scale at day 28 by treatment group, antiviral domain, FAS-ITT.
save_tex_table (
make_summary_table_antiviral (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-2-antiviral-summary" ))
make_summary_table_antiviral (fas_itt_dat)
```
```{r}
#| label: fig-who-antiviral-fas-itt
#| fig-cap: Day 28 WHO score, antiviral domain, FAS-ITT.
p1 <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
group_by (AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (AAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Antiviral intervention" )
p2 <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
AAssignment = factor (
AAssignment,
labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (AAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Antiviral intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-2-antiviral.pdf" ), p2, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-who-antiviral-avs-itt
#| fig-cap: Day 28 WHO score, antiviral domain, AVS-ITT.
lab <- c ("Usual care" , "Nafamostat" )
p1 <- avs_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
group_by (AAssignment = factor (AAssignment, labels = lab)) %>%
summarise (n = n ()) %>%
ggplot (., aes (AAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Antiviral intervention" )
p2 <- avs_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = lab),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (AAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Antiviral intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
fpth <- file.path (pth, "7-2-antiviral-avs-itt.pdf" )
ggsave (fpth, p2 + coord_flip (), height = 2.5 , width = 6 )
system (sprintf ("pdftoppm %s %s -png" , fpth, gsub (".pdf" , "" , fpth)))
p
```
```{r}
#| label: fig-who-antiviral-acs-itt
#| fig-cap: Day 28 WHO score, antiviral domain, ACS-ITT.
p2 <- acs_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
AAssignment = factor (
AAssignment,
labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (AAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-2-antiviral-acs-itt.pdf" ), p2, height = 3 , width = 6 )
p2
```
## Age
```{r}
#| label: fig-7-2-age
#| fig-cap: |
#| Distribution of WHO scale at day 28 by age group, FAS-ITT
#| fig-height: 4
pdat <- all_data %>%
filter_fas_itt () %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
who = ordered (D28_who, levels = 1 : 8 ),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = who)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Sex
```{r}
#| label: fig-7-2-sex
#| fig-cap: |
#| Distribution of WHO scale at day 28 by supplemental sex requirement, FAS-ITT
#| fig-height: 4
pdat <- all_data %>%
filter_fas_itt () %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
Sex,
who = ordered (D28_who, levels = 1 : 8 ),
.drop = F) %>%
group_by (Sex) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (Sex) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Sex, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
labs (y = "Number of \n participants" )
p2 <- ggplot (pdat, aes (Sex, p)) +
geom_col (aes (fill = who)) +
labs (y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-sex.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Oxygen
```{r}
#| label: fig-7-2-oxygen
#| fig-cap: |
#| Distribution of WHO scale at day 28 by supplemental oxygen requirement, FAS-ITT
#| fig-height: 4
pdat <- all_data %>%
filter_fas_itt () %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
supp_oxy = fct_explicit_na (factor (supp_oxy2, levels = 0 : 1 , labels = c ("Not required" , "Required" ))),
who = ordered (D28_who, levels = 1 : 8 ),
.drop = F) %>%
group_by (supp_oxy) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (supp_oxy) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (supp_oxy, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
labs (y = "Number of \n participants" ,
x = "Supplemental oxygen" )
p2 <- ggplot (pdat, aes (supp_oxy, p)) +
geom_col (aes (fill = who)) +
labs (x = "Supplemental oxygen" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-oxygen.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Country
```{r}
#| label: fig-7-2-country
#| fig-cap: |
#| Distribution of WHO scale at day 28 by country, FAS-ITT.
#| fig-height: 4
pdat <- all_data %>%
filter_fas_itt () %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
who = ordered (D28_who, levels = 1 : 8 ),
.drop = F) %>%
group_by (Country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (Country, p)) +
geom_col (aes (fill = who)) +
labs (x = "Country" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Site
```{r}
#| label: fig-7-2-site
#| fig-cap: |
#| Distribution of WHO scale by study site, FAS-ITT.
#| fig-height: 4
pdat <- all_data %>%
filter_fas_itt () %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (Location),
who = ordered (D28_who, levels = 1 : 8 )) %>%
complete (who = ordered (1 : 8 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = who)) +
labs (x = "Study site" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F, ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-7-2-calendar
#| fig-cap: |
#| Distribution of WHO scale by calendar time, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
who = ordered (D28_who, levels = 1 : 8 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = who)) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("WHO scale" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-2-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Sample Cumulative Logits
Proportional odds looks reasonable for all logits except 1.
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits, anticoagulation.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (CAssignment, D28_who) %>%
complete (CAssignment, D28_who, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (D28_who) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (D28_who != 28 )
ggplot (trt_logit, aes (D28_who, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits, antiviral.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (AAssignment, D28_who) %>%
complete (AAssignment, D28_who, fill = list (n = 0 )) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (AAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (D28_who) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (D28_who != 28 )
ggplot (trt_logit, aes (D28_who, rel_clogit)) +
facet_wrap ( ~ AAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
# Modelling {#modelling}
## FAS-ITT
```{r}
#| label: fit-model-fas-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
fas_itt_nona_dat,
ordmod,
outcome = "D28_who" ,
intercept = FALSE
)
# res2 <- fit_primary_model2(
# fas_itt_dat |> filter(!is.na(D28_who_low)),
# ordmodcens,
# outcome = c("D28_who_low", "D28_who_high"),
# intercept = FALSE
# )
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-fas-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), FAS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-2-primary-model-fas-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-fas-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-2-primary-model-fas-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio posterior summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-2-primary-model-epoch-site-terms-fas-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Posterior Predictive
```{r}
#| code-summary: Posterior predictive figure
#| label: fig-fas-itt-ppc
#| fig-cap: Posterior predictive check, FAS-ITT.
y_ppc <- res$ drws$ y_ppc
ppc_dat <- bind_cols (fas_itt_nona_dat, tibble (y_ppc = y_ppc))
grp_ppc1 <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
y_1 = mean (D28_who == 1 ),
ypp_1 = rvar_mean (y_ppc == 1 ),
y_2 = mean (D28_who <= 2 ),
ypp_2 = rvar_mean (y_ppc <= 2 ),
y_3 = mean (D28_who <= 3 ),
ypp_3 = rvar_mean (y_ppc <= 3 ),
y_4 = mean (D28_who <= 4 ),
ypp_4 = rvar_mean (y_ppc <= 4 ),
y_5 = mean (D28_who <= 5 ),
ypp_5 = rvar_mean (y_ppc <= 5 ),
y_6 = mean (D28_who <= 6 ),
ypp_6 = rvar_mean (y_ppc <= 6 ),
y_7 = mean (D28_who <= 7 ),
ypp_7 = rvar_mean (y_ppc <= 7 ),
y_8 = mean (D28_who <= 8 ),
ypp_8 = rvar_mean (y_ppc <= 8 )
) %>%
pivot_longer (y_1: ypp_8, names_to = c ("response" , "who" ), names_sep = "_" , values_to = "posterior" ) %>%
mutate (who = as.numeric (who))
}
grp_ppc2 <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
y_1 = mean (D28_who == 1 ),
ypp_1 = rvar_mean (y_ppc == 1 ),
y_2 = mean (D28_who == 2 ),
ypp_2 = rvar_mean (y_ppc == 2 ),
y_3 = mean (D28_who == 3 ),
ypp_3 = rvar_mean (y_ppc == 3 ),
y_4 = mean (D28_who == 4 ),
ypp_4 = rvar_mean (y_ppc == 4 ),
y_5 = mean (D28_who == 5 ),
ypp_5 = rvar_mean (y_ppc == 5 ),
y_6 = mean (D28_who == 6 ),
ypp_6 = rvar_mean (y_ppc == 6 ),
y_7 = mean (D28_who == 7 ),
ypp_7 = rvar_mean (y_ppc == 7 ),
y_8 = mean (D28_who == 8 ),
ypp_8 = rvar_mean (y_ppc == 8 )
) %>%
pivot_longer (y_1: ypp_8, names_to = c ("response" , "who" ),
names_sep = "_" , values_to = "posterior" ) %>%
mutate (who = as.numeric (who))
}
plot_grp_ppc <- function (dat) {
ggplot (dat %>% filter (response == "ypp" ), aes (x = who)) +
facet_wrap ( ~ grp, nrow = 1 ) +
stat_slabinterval (aes (ydist = posterior)) +
geom_point (data = dat %>% filter (response == "y" ),
aes (x = who, y = mean (posterior)),
colour = "red" ,
shape = 23 ) +
labs (x = "WHO outcome" ,
y = "Probability" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>% filter (response == "ypp" ), aes (y = who)) +
facet_wrap ( ~ grp, nrow = 1 ) +
stat_interval (aes (xdist = posterior), size = 2 ) +
geom_point (data = dat %>% filter (response == "y" ),
aes (y = who, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (xlab, breaks = c (0 , 0.5 ),
sec.axis = sec_axis (~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("WHO \n outcome" , breaks = c (1 ,3 ,5 ,7 )) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_epoch <- grp_ppc2 (sym ("epoch" )) %>%
mutate (grp = fct_inorder (factor (grp)))
pp_A <- grp_ppc2 (sym ("AAssignment" ))
pp_C <- grp_ppc2 (sym ("CAssignment" ))
pp_ctry <- grp_ppc2 (sym ("Country" ))
pp_site <- grp_ppc2 (sym ("site" )) %>%
left_join (ppc_dat %>% dplyr:: count (site, Country), by = c ("grp" = "site" ))
p0 <- plot_grp_ppc (pp_A, "Antiviral" , "" )
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (pp_site %>% filter (Country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (pp_site %>% filter (Country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (pp_site %>% filter (Country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (pp_site %>% filter (Country == "NZ" ), "Sites New Zealand" , "" )
p <- (p0 | p1 | p2) / p3 / p4 / p5 / (p6 | p7)+
plot_layout (guides = "collect" ) &
theme (legend.position = "bottom" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,
"7-2-primary-model-fas-itt-ppc.pdf" )
ggsave (pth, p, width = 6 , height = 5.75 )
p
```
## ACS-ITT
```{r}
#| label: fit-model-acs-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
acs_itt_nona_dat,
ordmod,
outcome = "D28_who" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-acs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), ACS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-2-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-acs-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts, ACS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-2-primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| label: fig-epoch-site-acs-itt-model
#| code-summary: Odds ratio posterior summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios, ACS-ITT.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-2-primary-model-epoch-site-terms-acs-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
## AVS-ITT
```{r}
p1 <- avs_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
group_by (CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (CAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Anticoagulation intervention" )
p2 <- fas_itt_nona_dat %>%
filter (! is.na (D28_who)) %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
D28_who = ordered (as.integer (D28_who), levels = 1 : 8 )
) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (CAssignment, p)) +
geom_col (aes (fill = D28_who)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "WHO scale" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-2-anticoagulation-avs-itt.pdf" ), p2, height = 3 , width = 6 )
p2
```
### Pre-specified model
- excludes epoch
```{r}
#| label: fit-model-avs-itt-pre-spec
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod_site,
outcome = "D28_who" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" , "ctry" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 1 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-pre-spec
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-2-primary-model-avs-itt-summary-table-pre-spec" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-pre-spec
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-2-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_site_terms (
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-2-primary-model-site-terms-avs-itt-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" , "gamma_site" , "tau_site" )) |> print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["tau_site" ])
```
:::
### Reduced Model
- excludes epoch, region, and site
```{r}
res <- ASCOTr:: fit_primary_model (
avs_itt_nona_dat,
ordmod0,
outcome = "D28_who" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 ),
intercept = FALSE
)
# res2 <- fit_primary_model2(
# avs_itt_dat |> filter(!is.na(D28_who_low)),
# ordmodcens0,
# outcome = c("D28_who_low", "D28_who_high"),
# vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"),
# beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5),
# intercept = FALSE
# )
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Oxygen requirement" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-2-primary-model-avs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-2-primary-model-avs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" ))
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
### Country
Rather than group Australia and New Zealand, keep them as separate countries in the model.
- excludes epoch
```{r}
#| label: fit-model-avs-itt-ctry2
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod_site,
outcome = "D28_who" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" , "ctry2" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 1 , 1 ),
intercept = FALSE ,
ctry_var = "ctry2_num" ,
site_var = "site2_num"
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" , "Nepal" , "New Zealand" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-ctry2
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-2-primary-model-avs-itt-summary-table-ctry2" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-ctry2
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-2-primary-model-avs-itt-odds-ratio-densities-ctry2.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_site_terms (
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("Australia" , "Nepal" , "New Zealand" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-2-primary-model-site-terms-avs-itt-ctry2.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" , "gamma_site" , "tau_site" )) |> print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["tau_site" ])
```
:::